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ABSTRACT 



Evolutions of X-ray clusters of galaxies are studied by N-body (shell model) 
+ mesh code (TVD) simulations on the assumption of spherical symmetry. We 
consider a density perturbation of 10^^ Mq composed of dark matter and gas in 
cold dark matter dominated universe with the cosmological density parameter, 
r2o = 1 or 0.2. A shock front appears during its initial collapse, moving outward 
as ambient gas accretes towards cluster center. The shock front separates the 
inner X-ray emitting, hot region, where gas is almost in hydrostatic equilibrium 
but with small radial infall (~ 100km s"^) being left, from the outer cool region, 
where gas falls almost freely and emits no X-rays. Gas inside the shock is 
strongly compressed and heated by shock so that X-ray luminosity rapidly rises 
in the early stage (until temperature reaches about virial). In the late stage, 
however, the X-ray luminosity rises only gradually due partly to the expansion of 
the inner high temperature region and partly to the increase of X-ray emissivity 
of gas as the results of continuous adiabatic compression inside the shock. We 
also find for clusters in lower density universe that the density distribution is 
generally less concentrated and, hence. X-ray luminosity more slowly rises than 
in higher density universe. 

The shock front structure, which was not clearly resolved in the previous 
SPH simulations, is clearly captured by the present simulations. Our results 
confirm that shock heating plays an important role in the heating process of 
intracluster medium. In addition, we find a sound wave propagating outward, 
thereby producing spatial modulations with amplitudes of ~ 10 % in the radial 
temperature and density profiles and time variations in the strength of the 
shock. Such modulations, if observed, could be used as a probe to investigate 
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the structure of clusters. 

Subject headings: galaxies: clustering — galaxies: intergalactic medium 
galaxies: X-ray — hydrodynamics 
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1. INTRODUCTION 



Clusters of galaxies (CGs) are luminous X-ray sources (Sarazin 1988). X-rays are 
thought to be emitted through thermal bremsstrahlung from very hot, optically thin plasma 
gas in intergalactic space, which is called intracluster medium (ICM). Recent X-ray and 
optical observations have revealed dynamical aspects of CGs (Fabricant, Kent & Kurtz 
1989; Kneib et al. 1995). Time evolution of the cluster X-ray luminosity function has been 
confirmed; each CG became brighter and the number density of CGs with a fixed X-ray 
luminosity increased through mergers towards zero redshift (Edge et al. 1990; Gioia et 
al.l990; Henry et al. 1992; Bower et al. 1993; Castander et al. 1994). To understand the 
physical structure of CGs, therefore, it is essential to consider non-equilibrium processes 
involved with the dynamical evolution of CGs. 

Recent extensive X-ray observations have established that the temperature profiles 
of ICM are nearly fiat in many clusters (Fabian 1994; Ohashi et al 1996). Although 
there is no widely accepted explanation for this, it is, at least, reasonable to expect that 
shock heating through gravitational collapse plays an important role. It then follows that 
the isothermality of gas may refiect the shape and time evolution of the gravitational 
potential well. What then governs the total dynamics of CGs? It is galaxies and dark 
matter (DM), both of which can be regarded as a coUisionless self-gravitating many-body 
system. Dynamics of ICM, collisional fiuid on a CG scale, simply follows the dynamical 
evolution of the gravitational potential exerted mainly by DM and galaxies. Thus, the 
thermal history of ICM is rather sensitive to the physical processes involving violent time 
variation of gravitational potential field (such as formation, merger, etc). Violent relaxation 
(Lynden-Bell 1967) is thought to play a crucial role there, although its physical significance 
is still poorly understood (Funato, Makino & Ebisuzaki 1992a,b; Takizawa & Inagaki 1996). 



-5- 

CGs were believed to be formed from overdense perturbations in the universe. They 
grew through gravitational instability and collapsed at z ~ 1. CGs have been growing 
by accreting ambient matter, which is still an ongoing process at the present. Gunn & 
Gott (1972) quantitatively discussed the growth of a spherical symmetric perturbation 
consisting of only collisionless particle in an expanding universe. Bertschinger (1985) found 
the self-similar solution describing evolution of spherical density fluctuation consisting of 
DM and gas in the Einstein-de Sitter universe (where the cosmological density parameter is 
fio = 1 and the cosmological constant is Aq = 0). 

There have been plenty of numerical studies performed so far regarding the formation 
and evolution of CGs by using N-body and hydrodynamic codes. Perrenod (1978) was 
the first to calculate the evolution of a spherical symmetric cluster with standard mesh 
hydrodynamic code. Three-dimensional calculations have been carried out recently, 
which mostly use smooth particle hydrodynamics (SPH) codes (Evrard 1990; Thomas & 
Couchman 1992; Bryan et al. 1994; Metzler & Evrard 1994; Navarro, Frenk & White 
1995, hereafter NEW). According to Evrard (1990), a shock front moved outward when 
gas around the cluster accreted onto the cluster center and a relatively fiat temperature 
profile was realized within the shock front. This result was significant in the sense that 
the above expectation was confirmed qualitatively. However, quantitative estimation may 
be problematic because of limited spatial resolutions and poor reliability of SPH code 
in calculations of shocks. It might be kept in mind that SPH codes are easy to work 
with and could give reasonable accuracy, however, they are not better for problems with 
discontinuities (such as strong shock and contact discontinuity) than mesh codes (Monaghan 
1992). On a larger scale, numerical simulations of cosmological structure formation by 
using N-body and hydrodynamic mesh code have been carried out mainly to investigate the 
statistical properties of CGs (Cen & Ostriker 1994; Kang et al. 1994; Anninos & Norman 
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1996). 

In this paper, we focus on the dynamical aspects of ICM, such as time- dependent 
properties of shock waves and the effects of shock heating on the evolution of CGs. For this 
purpose, we perform numerical simulations of spherical CGs with N-body (shell model) + 
mesh code (TVD). Note that the TVD code is one of the most useful tools to deal with 
spatial discontinuities. Since we assume spherical symmetry, the problem can be reduced to 
one-dimensional. Better spatial resolutions can be achieved, therefore. Using these codes, 
we calculate dynamical evolutions of a density perturbation of IO^^Mq consisting of DM 
and gas and collapsing at z ~ 1 in universe with fig = 1 or 0.2. It is also interesting to 
investigate how different cosmological models affect evolution of CGs. 

The rest of this paper is organized as follows. In section ^ we describe our numerical 
methods and the adopted initial conditions. In section |^ we present our results, discussing 
physical processes underlying the structural evolution of calculated CGs. In section ^ we 
summarize our results and discuss their implications. 

2. THE SIMULATIONS 

In the present study we regard CG consisting of two components: coUisionless particles 
corresponding to galaxies and DM, and gas corresponding to ICM. When calculating gravity 
both components are considered, although the former component dominates over the latter. 
We also assume spherical symmetry in all the calculations. 
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2.1. Basic Equations for Collisionless Particles 

For calculations of collisionless particles we adopt a shell model (Henon 1964). The 
distribution function of spherical symmetric stellar systems can be expressed as f{r,u,v,t) 
, where r is radial distance, u is radial velocity, v is tangential velocity, and t is time, 
respectively. So the state of the system is represented by N points in the (r, u, v) space 
(with N being the number of shells) and a trajectory of the i-th point at {ri,Ui,Vi) is 
calculated according to the equations of motion; 

(1) 
(2) 

for i = 1 r^ N, where G is the gravitational constant, Ai = ViVi (= constant in time) is 
angular momentum of the i-th shell, and M, is total mass (also including mass of gas) 
interior to Tj, respectively. Since it is convenient to carry out numerical calculations using 
comoving coordinates for our purpose, we transform (r^, Ui, Vi) to (i?j, Ui, Vi) as follows; 

Vi = a{t)Ri, (3) 

Ui = dRi + Ui, (4) 

where a{t) is the dimensionless scale factor of the universe and a represents the derivative 
of a{t) with respect to time. Equations (P and (|^) are then transformed into 

oiXj Ui , , 

dt ~ {aR,y {aRiY a' '' ^ ' 

respectively, by using equations (|]) and (^). Equations (H) and (||) are integrated by 
using leap-frog method. As to the inner boundary condition we set a reflecting wall at 
Trw = 0.02/(1 + z) and impose that when a shell reaches the wall, r, < r™, that shell is 



elastically reflected. Only shells with rather small angular momentum is influenced by the 
wall. 

2.2. Basic Equations for Gas 

For gas components basic equations in the comoving frame are as follows, 
dp 1 d a 2 






^ ^^■p\ , ^ ^ ^„zj„. \ _ '^/^^™,2 37 ^\ ^pHvg^s 



where p is gas density, Wgas is radial velocity of gas, P is gas pressure, and 7 is the adiabatic 
exponent, respectively. The total energy of gas per unit mass, E, and the enthalpy per unit 
mass, H, are given by 

^ - M^^f ■ ('«' 

H = ^_- + ^^, (11) 

respectively, and qr is deflned by 

^^--(^-ai?, (12) 

where, M^j is the total mass (including gas and colhsionless particles) inside R. 

We neglect viscosity and angular momentum of gas. We also assume that gas is ideal 
gas with 7 = 5/3. A second order up-wind TVD code (minmod limiter) is used for our 
simulations (Hirsch 1990). Number of mesh point is 500 and one mesh spacing corresponds 
to Ar = 0.02/(1 + z) Mpc with z being a redshift [1 + z = a(0)/a(t)]. We set boundary 
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condition as follows. The inner edge is assumed to be a perfectly reflecting point; 

P-i = Po = Pi, (13) 

P_i = Po = Pi, (14) 

^^gas -1 = ^^gas,0 = 0, (15) 

where pi, Pi and fgas,i is gas density, gas pressure and radial velocity of gas at the i-th mesh 
point, and the first mesh point corresponds to the inner boundary. The 0-th and —1-st 
points are necessary for calculations (to derive spatial derivatives of physical quantities). 
The outer edge is assumed to be a perfectly transmitting surface; 

g„+2 = 3g„ - 2g„_i, (16) 

g„+i = 2g„ - g„_i, (17) 

where g^ is any physical quantity of gas at the i-th mesh point and the n-th corresponds to 
the outer boundary. Again, the (n+ l)-th and (n + 2)-th points are necessary for calculation 
purpose. 

2.3. Models and Initial Conditions 

In this paper, calculations are carried out from Zi^i = 10 to the present time {zq = 0). 
We adopt cosmological models with no cosmological constant, Aq = 0. When fig = 1, 
therefore, we find 

t \l 



3^0 



«(^)=(t77^ ' (1^ 



while when < r2n < 1 we have 



o(«) = ^,^"° (cosh-tf.-!). (19) 
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where Hq is the Hubble constant and is set to be Hq = 100 (km s~^ Mpc^^) in 
transformations of length and time. 

We make initial density profiles in the same manner as Peebles (1982). At first we 
prepare A^ concentric shells with uniform density being equal to the mean density of the 
universe at Zi^^i = 10. Then a density fiuctuation is introduced by perturbing the radius and 
velocity of each shell as 



(0) 

Ti = r) 



.(0)^ 



l + Krr) , (21) 



where r] is the unperturbed coordinate and 5{r) represents the perturbation as a function 
of r (specified below). The velocity perturbation is written, using Zel'dovich approximation, 
by 



Ui = H{zi^i)rf^ 



.(0)^ 



l + 26(rn , (22) 



where H{zi^i) is the Hubble constant at the initial epoch, and we used 5 oc t^^^; the 
relation which holds exactly in the Einstein-de Sitter universe. The functional form of the 
perturbation is assumed as 

f-%cos2(^^) (0<r<ro) 
[ (ro < r) 

where tq and So are, respectively, the initial size of the fiuctuation on the comoving scale 
and the parameter for displacement (see Nakamura 1996). We take {SqjTq) = (0.4,9.0 
Mpc) for Model A (with fio = 1-0) and (5o,ro) = (0.5,15.4 Mpc) for Model B (with 
Qq = 0.2), respectively. In both models the fiuctuation contains a mass of about 10^^ Mq 
and correspond to typical density peaks with ~ 1.5o" in the CDM power spectrum with 
(Tg = 0.96. Here, a is the rms density fiuctuations on a cluster scale and ag is the rms of 
mass fiuctuation on scale 8 Mpc; ag = 0.96 was obtained from the observation of nearby 
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galaxy distribution (see Suto 1993). The profile of tlie initial density perturbation and the 
ratio of ((5M/M) /erg cdm are illustrated in figure ^for Model A (by the solid line) and B (by 
the dotted line), respectively. Here, {6M/M) is the mass fluctuation averaged over a scale r 
,whereas cr8,cDM is the rms of mass fluctuation averaged over a scale r obtained from CDM 
power spectrum with as = 0.96. 

Fig. 1.- 

The initial conditions of gas are set as follows. At first density of gas was everywhere 
set to be 1/10 of the mean density of the universe at z = 10 and temperature of gas (Tgas,i) 
was everywhere taken to be 10''' K except in Model LT, where we assumed Tgas,i = lO^K 
initially. Then adiabatic fluctuation (cf. Eq. 16 and 17) was imposed so that gas density 
became always 1/10 of that of DM. Note that gas temperature distribution becomes 
nonuniform, accordingly. We calculated 3 models in total. Model parameters are listed in 
table [I|. 

3. RESULTS 
3.1. Overview of Evolution 

At first we overview the evolution of Model A cluster. The evolution of density 
and radial velocity profiles of DM are displayed in figure 0. The different types of lines 
correspond to different redshifts (and different times): z = 2.5 (t = 1 Gyr) by the long dash 
line; z = 1.2 {t = 2 Gyr) by the short dash line; z = 0.4 (t = 4 Gyr) by the dotted line; and 
z = (t = 6.5 Gyr) by the solid line, respectively. 

Fig. 2.— 
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Density profile basically obeys a power law, p^M oc r~^'^, and evolves in a self-similar 
fashion after z ~ 1. The evolutions in the radial distributions of radial velocity, density, 
temperature, pressure, and entropy (each of the gas component) are summarized in figure 
^. The different types of lines correspond to same redshifts as in figure 0. 

Fig. 3.- 

There are fundamental features commonly seen in the course of evolution of all the 
models, which can be summarized as follows. 

1. Before the initiation of DM collapse (at z > 2.5; t < 1 Gyr) gas continues to expand, 
following the cosmological expansion. When DM begins to collapse, a shock wave 
emerges in the central part and moves outwards, accreting ambient gas towards the 
center. 

2. The shock front separates the inner hot region from the outer cool region. In the 
inner region, gas is almost in hydrostatic equilibrium, although bulk velocity of radial 
infall~ 100 km s~^ still remains. The temperature profile is relatively fiat and gas is 
hot enough to emit X-ray. In the outer region, in contrast, gas falls almost freely and 
is too cool to emit X-ray. 

3. Density profile of gas evolves self-similarly as ngas oc r~^'^ inside the shock front except 
near the center, where density profile is rather fiat. Even after the passage of the 
shock, density, temperature, and, therefore, pressure, gradually increase with time. 
This indicates that the inner region is not perfectly in hydrostatic equilibrium (will 
be discussed in §3.3). 

4. Entropy profile shows an overall increase outwards, suggesting larger entropy 
production taking place as the shock moves outward. There are wavy features seen. 
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especially, in the distributions of temperature and entropy. These seem to be related 
to sound wave propagation (will be discussed in §3.4). 

Let us next examine each item in more details and discuss similarities and differences 
between different models. 

3.2. Density Profiles 

Density profiles of DM aX z = can well be fitted with the /?-model; 

2"! — 3/3dm/2 



Pdm(^) = Pdm.o 






(24) 



Here, tdm and /?dm are fitting parameters. We fit the results of simulations inside rioo, the 
radius where the mean interior density is 100 times of critical density at z = 0. The results 
of the fitting are summarized in table |^. Note that these results may depend on initial 
density profile. 

In the same way density profiles of gas at 2; = inside the shock fronts can be fitted 
with the /3-model, 



'^gasv j ^gas,0 



^<H 



f \ 2"| 'jPgas/'^ 



^gas 



(25) 



where rgas and /3gas are fitting parameters. The results are listed in table 0. In all models 
we find /5dm ~ Pgas ~ 0.9 and tdm ~ ^gas- 

We also depict nondimensional density profiles at 2; = in figure ^ for Model A (by the 
solid line) and Model B (by the dotted line), respectively, where density is scaled with pco, 
critical density of the universe at 2; = 0, and radius is scaled with rioo. Nondimensional DM 
density profiles look similar among two models, whereas gas component expands slightly in 
Model B, compared with that in Model A. 
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Fig. 4.— 

Let us finally compare our results with the previous results (NFW). As for the density 
profiles of the clusters in Einstein-de Sitter Universe (Models A and LT), the core radii 
are smaller than those of NFW (who obtained tdm ~ '"gas ~ 0.2), while /5dm and /3gas are 
similar. On the other hand, central DM density, Pdm,0; is similar to that of NFW, whereas 
central gas density, Pgas,o, of model A is lower than that of NFW (our model LT is consistent 
with their model; ngas,o ~ 5 x 10^^). This can be explained in terms of different initial gas 



temperatures (see section p.7|) . 



3.3. Temperature Profiles 

Temperature and entropy profiles are shown in figure |^ in the upper and lower panels, 
respectively. The different types of lines correspond to different models: Model A by the 
solid hne; Model B by the dotted line; and Model LT by the short dash line, respectively. 
There are common characters as follows. Temperatures gradually fall outwards until the 
shock front. There are small temperature fluctuations seen. If we evaluate the propagation 
speed of fluctuation pattern, it is ~ 700km/s, of the order of the sound velocity inside 
the shock. In addition, entropy pattern does not change much during the course of wave 
propagation, indicating that structural variation is adiabatic. We may thus conclude that 
the temperature fluctuation arises due to sound wave propagation. 

Inside the shock fronts temperature is nearly virial, high enough to emit X-ray. Entropy 
monotonically rises from the center to the shock front. That is, near the center gas is 
heated up mainly through adiabatic compression, thus possessing relatively lower entropy. 
Gas in the outer region, on the other hand, will eventually be heated through shocks, which 
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effectively transform kinetic energy of accreting gas into thermal energy at their surface. 
This explains relatively large entropy just inside the shock. Note that regions with entropy 
decreasing outward are convectively unstable. Convective motions, if occur, will smear out 
such a feature. 

Fig. 5.— 

To understand why temperature, density and entropy steadily increase with time even 
after the passage of the shock front and how sound waves arise and propagate outward, we 
check the balance between pressure gradient and gravity in the inner parts. Figure |^ shows 
the ratio of the absolute value of pressure gradient force, \dP/dr\, and gravitational force 
inside the shock surface, Fg = GM^pgas/T^, from t = 5.5 Gyr to 6.5 Gyr. If the system is 
perfectly in hydrostatic equilibrium, \dP/dr\/Fg should be unity. 

This figure shows a clear tendency that gravitational force overcomes pressure gradient 
force inside the shock front, thus inducing radial gas inflow from outside. Therefore, gas 
is adiabatically compressed by the infalling material from outside. Importantly, the ratio 
changes with time. Since DM density profile hardly changes with time near the core, so 
does the gravity force; the ratio changes are purely due to time and spatial variations in 
pressure profiles. When ambient gas suddenly falls towards the center, pressure at the 
core abruptly increases. This gives rise to an outwardly propagating sound wave, since the 
central point is a reflecting boundary due to spherial symmetry (cf. §2.2). Figure |^ shows 
radial velocity profiles of gas inside the shock from t = 5.5 Gyr to 6.5 Gyr. We confirm that 
radial infall systematically remains even after the passage of shock. In addition we confirm 
fluctuation pattern propagating outwards with a speed roughly equal to the sound velocity. 

Fig. 6.— 
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Fig. 7.- 

In our simulations thermal conduction is neglected. It is possible, however, that 
thermal conduction, if efficient, will erase temperature fluctuations on a small scale. It is 
thus worthwhile evaluating conduction timescale in the simulated clusters. Conduction 
time scale is generally expressed as (cf. Sarazin 1988) 

tcond ~ — , 26 

K 

where rie is electron number density, /t is scale length of temperature gradient, k-Q is the 
Boltzmann constant, and thermal conductivity for hydrogen plasma is (Spitzer 1962) 

'^ ^ 4.6 X 10^3 (^^j (it) (^^g^"'^^''K-i), (27) 

where In A, Coulomb logarithm, is 



lnA = 37.8 + ln 



From equations (^) and (p7|) , we derive 



108KAio-3cm-3 



/ Te V^^V ^T \VlnA 



(2^ 



'— -K,5.^)(^)-'(e|;,)(^). m 

For Model A, for example, we find t^ond ~ 1-5 x lO^yr at r ^ 0.1 Mpc (where 
rie ~ IQ-^cm-^, Te ~ 5 X lO'^ K, and /t ~ 0.1 Mpc) and tcond ~ 1-5 x lOV at r ^ IMpc 
(where rie ~ 10~^cm"'^, To ~ 2 x 10^ K, and Zt ~ 0.1 Mpc). At both radii, hence, conduction 
seems to be efficient. The same is true for Model B. However, we should note that the 
usage of the classical conductivity (Spitzer 1962) is in question for CGs. In fact, we cannot 
explain the existence of cooling flows as long as we employ the classical one (Binney & 
Cowie 1981). Instead, it is suggested that tangled magnetic field (Rosner & Tucker 1989) 
or plasma instabilities (Pistinner & Shaviv 1996) are likely to suppress heat conduction 
significantly in CGs. Temperature fluctuations can then survive. 
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3.4. Evolution of Shock Surface 

As we have seen in figure |^, sound wave propagates outward, modulating temperature 
profile. This also affects shock front propagation, since sound velocity inside the shock is 
three times greater than the front velocity in strong shock limits. We illustrate in figure ^ 
the radius of a shock surface (rghock) as a function of the look-back time, time measured 
from the present time; i.e., tib = corresponds to the present time {z = 0). The different 
types of lines correspond to same models as in figure |^. Shock surface moves with nearly 
a constant velocity (~ 200km s^^). The radius of shock surface at z = 0, r^, and mean 
propagation velocity, Va, are listed in table |^. Again, there are wavy features seen in this 
figure. 

Fig. 8.— 

To understand the physics causing modulating features, we plot time variation of shock 
strength and shock radius (rshock) for Model A in figure |^. To evaluate shock strength 
we use Av = vin — fout with Vin and fout being radial gas velocities inside and outside the 
shock surface, respectively, in the upper panel, and TjAS where Tj is the pre-shock gas 
temperature in the unit of lO'^K and AS = Sin — Sont with Sin and S'out being specific 
entropies inside and outside the shock surface, respectively, in the lower panel. Since this 
quantity is proportional to heat produced through shock heating, it is a good representation 
of the shock strength. Both panels show time modulation in shock strengths. This in turn 
creates spatial modulation in entropy profile, since shock radius moves outward with time; 
namely, at the radius over which the shock passed with its maximum (minimum) strength 
in time, radial entropy profile exhibits a rapid (or slow) rise outward. The time modulation 
in shock strengths is likely to be caused by shock radius oscillation, and this oscillation is, 
as we discussed above, caused by the sound wave propagation. Note that while the wavy 
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pattern in temperature profile at a fixed radius varies with time because of sound wave 
propagation, that in entropy profile does not, since entropy profile is unaffected by sound 
waves. 

Fiff. 9.— 



3.5. Evolution of X-Ray Luminosity 

Time evolutions of X-ray luminosity, L(tib), and the normalized luminosity, 
Ln(tib) = L{tYb) / L{tvo = 0), are plotted for each model in figure |l^ (a) and (b), respectively. 
When calculating luminosity, we assume thermal bremsstrahlung of optically thin plasma 
(Rybicki & Lightman 1979), 

^ff ^ dW ^ ^_^ ^ 10-27^1/2^2-^ (erg g-^cm-^), (30) 

dV at 

where T is gas temperature, n is number density, qb is a frequency average of velocity- 
averaged Gaunt factor. In this paper we set qb = 1-2. Emission from cool gas with 
temperature T < 10^ K was neglected because we are interested in time variation of X-ray 
luminosity. 

Fig. 10.— 

Figure |Tl| plot evolution of L and L^ against z instead of time. 
Fig. 11.— 

In all models the luminosity rapidly rises just before the appearance of a shock wave, 
and then rises gradually afterwards. This behavior is due to the two distinct phases of 
evolution of gas. 
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The first phase (before the emergence of shock wave): Through accretion of 

ambient gas towards the center gas in the central region is compressed adiabatically 
so that temperature and density there rapidly rise until temperature reaches about 
virial temperature. The rapid increase in L is thus due to a rapid rise in temperature 
of the central region. 

The second phase (after the appearance of shock wave): As the shock wave 

propagates outwards, high temperature region expands, thus increasing L. At the 
same time, gas inside the shock wave continues to be compressed adiabatically, 
thereby its emissivity being increased gradually. Luminosity thus increases faster than 

3.6. Dependence on the Cosmological Density Parameter 

In this subsection we compare the results of Model A {VLq = 1) and Model B {Qq = 0.2) 
to discuss the Qq dependence of the cluster evolution under the condition that the 
perturbation amplitudes at 2; = 10 and total masses contained in the perturbations are 
similar. 

From table |^, we see that central DM density, pDM,0)is proportional to Qq. The DM 
distribution is less concentrated in Model B than in Model A, because the ambient gas 
density is lower in a lower density universe so that a larger volume is needed to contain the 
same amount of mass (~ IO^^Mq) initially. The absolute values of DM density are thus 
different among these models, but the shape of the DM density profiles are similar. This 
explains why nondimensional DM density profiles look similar among both models (figure 
^). On the other hand, central density of gas, ngas,0) depends also on its initial entropy and 
is, hence, not strictly proportional to Qq (see table ^. This is responsible for the differences 
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in nondimensional gas density profiles among two models (figure ^). 

In both models X-ray luminosity rises, however, the increase in L{t) is slower in Model 
B than in Model A. This tendency is more clearly seen when L{t) is plotted against times, 
rather than against redshifts (see Fig. p!0|b, |ll|b). In the latter figure, the difference can be 
recognized at z ~ 0.5, but not at 2; > 1. This is because smaller amount of DM and gas 
surrounding a central condensation in a lower Qq universe. 

3.7. Dependence on the Initial Temperature of Gas 

We calculated a model with lower initial temperature (LT) to see how initial 
temperature affects the later evolution. The temperature and entropy profiles at z = of 
Models A and LT are depicted in figure |. There is no significant difference among these 
models, especially in the structure inside the shock front. The thermal property of gas 
outside the shock front has little influence on the thermal property of gas inside because of 
enormous entropy production at the shock surface. 

On the other hand, the central densities are different among these models (figure ^). 
According to rigas.o estimated from the fitting data of table |^ and figure ^, central density 
of Model LT is about three times greater than that of Model A. In the density profiles at 
r > 0.3 Mpc, however, any difference can hardly be seen. 

Fig. 12.— 

In the central region, gas is adiabatically compressed until the temperature reaches 
about the virial temperature. So lower initial gas temperature results in higher gas density, 
if the virial temperatures are the same in both models. However, the difference in central 
gas densities cannot be perfectly explained by this picture. If there is no entropy production 
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(i.e., gas is perfectly adiabatically compressed), central density of Model LT should be 
about thirty times greater than that of Model A. This means, shock heating took away the 
information about the initial gas temperature also in the central region. 

Because of higher central gas density Model LT cluster is more luminous in X-ray 
than Model A cluster (figure p!0|a, [Tl|b). Substantial difference, however, cannot be seen in 



evolutionary behavior of X-ray luminosity of both model, (figure p!0|b, |TT| 



4. SUMMARY AND DISCUSSION 

We carry out the numerical simulations of spherical clusters of galaxies with shell 
model for DM and second order up-wind TVD scheme for ICM to examine structural 
evolution of ICM. Shock front moves outwards as gas accretes towards cluster center, 
yielding a relatively fiat temperature profile inside the shock front. Density and pressure 
profiles evolves in a self-similar fashion. 

X-ray luminosity increases with time in two steps. At the initial collapse of DM gas 
in the central part is at first adiabatically compressed through accretion of ambient gas 
towards the center. Eventually, a shock wave appears near the center and X-ray luminosity 
rapidly rises until temperature increases and reaches about virial temperature via shock 
heating. In the late stage, in contrast, the luminosity rises only gradually, since the inner 
region already emits strong X-rays. The gradual brightening is due partly to the expansion 
of the high temperature region and partly to increasing X-ray emissivity of gas as the 
results of continuous adiabatic compression the gas inside the shock. 

If we compare two clusters with the same density fluctuation amplitudes at 2; = 10 and 
with the same total masses but in different mean-density unvierses, the DM distribution is 
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less concentrated in clusters in lower density universe. Hence, X-ray luminosity of clusters 
rises more slowly than in higher density universe. 

Initial gas temperature has some influence on the central gas density. Higher initial 
temperature results in lower central density. This corresponds to the case that reheating 
of the ICM by, e.g., proto-galaxies is substantial. Density distribution in the outer region 
and temperature proflle inside the shock front are, however, hardly influenced by changes 
in the initial temperature. Thus, the inclusion of reheating process modifles the scaling law 
between X-ray luminosity and temperature in the favorite way to reproduce the observed 
relation (NFW). Note that the epoch of the shock emergence depends on the initial speciflc 
entropy at the core; in the presence of reheating process (so that the initial gas temperature 
is ~ lO^K), speciflc entropy at the core is already high enough, and so the appearance of 
shock may be delayed. 

In table § we compare the time dependent properties of Model A cluster and those of 
the self-similar solution by Bertschinger (1985). Both behavior looks very similar. Note 
that the self-similarity can be seen in all the calculated models. Although Einstein-de Sitter 
universe is assumed in Bertschinger (1985), we flnd that the self-similarity can be found for 
cluster evolution in low density universe {Qq = 0.2, Aq = 0). 

Our results regarding the present proflles of density, temperature, and so on roughly 
coincide with those of the previous SPH simulations (Evrard 1990; NFW). However, the 
structure of the shock front, which was not well resolved in the previous SPH simulations, 
is now clearly captured in the present mesh-code simulations; our results show that hot 
gas of X-ray CG is separated with a deflnite boundary and that shock heating plays an 
important role for the heating process of ICM. We conflrmed the persistence of radial 
infall of gas after the passage of the shock front. This results in adiabatic compression of 
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the inner parts, inducing gradual temperature increase. We also found time variations in 
the strengths of the shock due to a sound wave propagation over the entire cluster, which 
modulates the radial distributions of temperature and density with relative amplitudes of 
about 10 %. Note that the property of the sound wave may depend on our assumption 
of spherical symmetry. It is open to question how such sound wave behaves in a realistic 
three-dimensional situation; i.e, when mass accretion takes place in a nonaxisymmetric way 
(e.g., by merger of multiple density condensations). If the temperature modulation will be 
observed with future X-ray mission (such as ASTRO-E), this can be used as a good probe 
to investigate the structure, especially, the mass distribution of CGs, since sound wave 
properties sensitively depend on the shape and the depth of the gravitational potential well, 
and thus on the DM mass distribution. 

We would like to thank F.E.Nakamura and M.Hattori for helpful discussions and 
comments. This work is in part supported by Research Fellowships of the Japan Society 
for the Promotion of Science for Young Scientists (M.T.), and by the Grants-in Aid of the 
Ministry of Education, Science, Sports, and Culture of Japan, 06233101, 08640329 (S.M.). 
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Model 


^^o 


^0 


(Mpc) 


Tgas,i (K) 


no 


A 


0.4 




9.0 


10^ 


1.0 


B 


0.5 




15.4 


10^ 


0.2 


LT 


0.4 




9.0 


10^ 


1.0 



Table 1: Parameters of each model. 
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Model 


Pdm,o (10-2' g 


cm 


-') 


roM (Mpc) 


/^DM 


^100 


A 


78.5 






0.067 


0.87 


1.35 


B 


14.9 






0.13 


0.93 


1.04 


LT 


83.3 






0.070 


0.87 


1.34 



Table 2: Density profiles of DM aX z = 0. 
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Model 


^gas,0 (10"^ 


cm~ 


-') 


rgas (Mpc) 


Pgas 


A 


9.9 






0.082 


0.80 


B 


3.9 






0.13 


0.93 


LT 


30 






0.065 


0.87 



Table 3: Density profiles of gas at 2; = 0. 
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Model 


Ta (Mpc) 


Va 


(km s-i) 


A 


1.2 




210 


B 


1.4 




230 


LT 


1.2 




230 



Table 4: The arrival radius of shock surface at z = 0, r^, and mean propagating velocity, Va 
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Pdm(^) Pgas(^) P{r) rshock(i) 



^ ^-2.6 ^-2.4 ^-3.0 ^0.83 

Bertschinger(1985) r-2-25 ^-2.4 ^-2.9 ^o.89 

Table 5: The comparison between the behavior of Al and the self-similar solution of 
Bertschinger (1985) 
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Figure Captions 

figure |l] (a) Radial profiles of initial density perturbations and (b) ((5M/M) /crg^cDM, where 
{6M/M) is the mass fluctuation averaged over a scale r and crg^cDM is the rms of 
mass fluctuation of scale r obtained from CDM power spectrum with cig being a 
normalization. The solid line correspond to Model A and the dotted line corresponds 
model B. 

figure ^ The evolution of (a) density profile and (b) radial velocity profile of DM of Model 
A cluster. The different lines correspond to different redshifts (and different times): 
z = 2.5 (t = 1 Gyr) by the long dash line; z = 1.2 (t = 2 Gyr) by the short dash line; 
z = 0.4 (t = 4 Gyr) by the dotted line; and z = {t = 6.5 Gyr) by the solid line, 
respectively. 

figure ^ Evolution of Model A cluster. Radial profiles of representative physical quantities 
of gas are depicted: (a) radial velocity, (b) density, (c) temperature, (d) pressure, and 
(e) entropy. The different types of lines correspond to same redshifts as in figure ^ 

figure ^ Nondimensional density profiles at z = for Models A (by the solid line) and for 
B (by the dotted line), respectively. Here, density is scaled with pco, critical density 
of the universe at z = 0, and radius is scaled with rioo- 

figure ^ (a) Temperature profiles and (b) entropy profiles of each model at 2; = 0. Entropy 
is normalized to zero at the outer boundary. The different types of lines correspond 
to different models: Model A by the solid line; Model B by the dotted line; and Model 
LT by the short dash line, respectively. 

figure ^ The ratio of pressure gradient force and gravitational force inside the shock 

surface of Model A cluster from t = 5.5 Gyr to 6.5 Gyr. If the system is perfectly in 
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hydrostatic equilibrium, the ratio, \dP/dr\/Fg, should be unity. 

figure ^ Radial velocity of gas inside shock front from t = 5.5 Gyr to 6.5 Gyr of model A. 
Sound wave propagates outwardly. Radial infall systematically remains even after the 
passage of the shock. 

figure ^ Time evolution of the radius of a shock surface, rshock- The abscissa is the 

look-back time, time measured from the present time; tib = at 2; = 0. Shock surface 
moves with a nearly constant velocity (~ 200 km s~^). The different types of lines 
correspond to same models as in figure |^. 

figure ^ Evolution of shock strength and the shock radius of Model A. (a) Time evolution 
of '"shock by the dotted line and that of Av, jump in the radial velocity of gas over the 
shock surface, (b) That of TiAS, where Tj is the pre-shock gas temperature in the 
unit of lO^K and AS* is the jump in the entropy over the shock surface. This quantity 
is proportional to heat produced through shock heating. 

figure |1^ (a) Time evolution of X-ray luminosity, L(tib), and (b) that of normalized 
luminosity, Ln(tib) = L(tib)/-L(tib = 0). In all models L rapidly rises just before the 
generation of a shock wave, and then rises gradually afterwards. The different types 
of lines correspond to same models as in figure ^. 



figure 11 Evolution of L and L^ plotted against z. 



figure |T^ Gas density profiles for Model A (solid line) and for Model LT (dotted line). 
Central density is about three times greater in Model LT cluster than in Model A 
cluster. In the region of r > 0.3Mpc, however, this difference can hardly be seen. 
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